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Abstract 

Our article considers the class of recently developed stochastic models that combine claims payments and incurred 
losses information into a coherent reserving methodology. In particular, we develop a family of hierarchical Bayesian 



12| and Gogol 



In the process 



paid-incurred-claims models, combining the claims reserving models of Hertig 

we extend the independent log-normal model of Merz and Wiithrich [l^ by incorporating different dependence 
structures using a Data- Augmented mixture Copula paid-incurred claims model. 

The usefulness of incorporating both payment and incurred losses into estimating of the full predictive distribution 
of the outstanding loss liabilities and the resulting reserves is demonstrated in the following cases: (i) an independent 
payment data model; (ii) the independent payment and incurred claims data model of Merz and Wiithrich [igI ]: 
(iii) a novel dependent lag- year telescoping block diagonal Gaussian copula payment and incurred claim data model 
incorporating conjugacy via transformation; (iv) a novel data-augmented mixture Archimedean copula dependent 
payment and incurred claim data model. 

Inference in such models is developed by adaptive Markov chain Monte Carlo sampling algorithms. These incor- 
porate a data-augmentation framework utilised to efficiently evaluate the likelihood for the copula based payment 
and incurred claim model in the loss reserving triangles. The adaptation strategy of the Markov chain Monte Carlo 
is based on two components. The first component uses an adaptive strategy for learning the posterior structures 
for the parameters defined over a Euclidean space and the second component deals with an adaptive learning of 
the posterior for the covariance matrices restricted to the Riemann manifold corresponding to the space of positive 
definite matrices for the linear dependence structure specified for the payment and incurred claim model. 

Key words: Chain Ladder Models, Claims Reserving, Data Augmentation, Adaptive Markov Chain Monte Carlo 



1. Introduction 

As discussed in Merz and Wiithrich [l^ the main task of reserving actuaries is to predict ultimate loss ratios 
and outstanding loss liabilities. In general such predictions are based on past information that comes from a 
variety of sources. Under a credibility based framework, the weighting of different data sources and their relative 
contribution to the estimated reserve is difficult to determine. Therefore, it is important to consider developing 
a unified prediction framework for the outstanding loss liabilities, known as the paid-incurred-claims (PIC) class 
of models. However, to date only simple dependence structures have been considered, with three parameters for 
the correlations which were not incorporated into the formal Bayesian estimation approach, and instead fixed 
deterministically a priori. There are two technical difficulties in extending the current restrictive assumptions 
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within a Bayesian framework. The first is being able to generate the positive definite matrices; the second is 
evakiating the joint likelihood of the mixture copula defined over the observed payments and incurred losses in 
each accident year row of the reserving matrix. Our article significantly extends the dependence structure of 
current PIC models by solving these two problems. The first problem is resolved through utilisation of a class of 
matrix-variate Inverse- Wishart priors coupled with an adaptive Markov chain sampler that restricts the proposed 
Markov chain states to remain on the manifold of such matrices. The second problem is solved by using a data 
augmentation strategy which treats the unobserved parts of the loss triangle as missing data so that one can 
perform evaluation of the copula based likelihood required for inference on the model parameters. 

In order to ensure the financial security of an insurance company, it is important to predict future claims liabilities 
and obtain the corresponding prediction intervals which take into account parameter uncertainty. The PIC model 
is a claims reserving method which statistically combines information about claims payments and incurred losses. 
It allows actuaries to best utilise the available information for loss reserves. The Munich chain ladder method 
introduced by Quarg and Mack 



24| is one of the first claims reserving approaches in the actuarial literature to 
unify outstanding loss liability prediction based on both sources of information. This method aims to reduce 
the gap between the two chain ladder predictions that are based on claims payments and incurred losses data, 
respectively. It is achieved by adjusting the chain ladder factors with paid-incurred ratios to reduce the gap 
between the two predictions. The main drawback with the Munich chain ladder method is that it involves several 
parameter estimates whose precisions are difficult to quantify within a stochastic model framework. 

Merz and Wiithrich recently introduced a log-normal PIC chain model and used Bayesian methods to estimate 
the missing (future) part of the claims reserving triangles based on both payment and loss incurred information. 
Its major advantage is that the full predictive distribution of the outstanding loss liabilities can be quantified. 
One important limitation of the model of Merz and Wiithrich [l^ is that it does not develop the dependence 
properties of the PIC model that will be applicable to loss reserving data observed in practice. Our article extends 
the proposed Bayesian PIC models of Merz and Wiithrich [l^ to capture additional dependence structures. 

1.1. Brief Background 

Dependence within payment data, within incurred loss data, and even between payment and incurred loss data 
commonly exists due to the nature of the loss process. Payment and incurred loss ratios in the previous development 
period are likely to impact that of the next development period. Hence, correlation between development periods 
is practically appealing in claims reserving practice. Moreover, incurred loss is essentially payment data plus 
case estimates which are projections foreseen by case managers to estimate the remaining payments. Correlation 
between payment and incurred loss data is also found. Happ and Wuthrich 



11 



propose a fixed covariance structure 



to describe the correlation between payment and incurred loss, assuming that the correlations between different 
development periods are identical. In reality, correlations differ across development periods for various reasons, such 
as different stages of the life cycle for a claim and internal policy changes. In order to fully incorporate the actual 
correlations, we introduce a block covariance structure to allow for the variation between different development 
periods within payment and incurred losses. We also develop a second class of hierarchical mixture of copulas 
models. 

The estimation challenge involves constructing and sampling from the resulting Bayesian models for PIC with 
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flexible dependence structures. To specify the model, we vectorize the triangular random structures for payments 
and incurred loss and, applying appropriate permutations, we then assume a copula dependence structure on the 
vectorized data. We use a Gaussian copula with an unknown correlation matrix, which is restricted to be block 
diagonal for parsimony, or a mixture Archimcdian copulas across development periods. 

We estimate the Bayesian models by MCMC methods, using data augmentation to generate missing data values 
in the loss triangle and use an adaptive Metropolis algorithm to generate the unknown parameters. Bayesian 
simulation methodology is used to carry out inference on all aspects of the models considered and to obtain 
predictive distributions for reserves. 



1.2. Contributions 

We design a novel class of PIC models and illustrate it with two examples. The first involves a mixture of Clayton 
and Gumbel copulas for upper and lower tail dependence features in the development years for payments and 
incurred losses. The second example involves a Gaussian copula model in which the covariance structure is a 
telescoping block diagonal form representation which captures dependence between development lag years in the 
payments and incurred losses. By a telescoping block diagonal matrix we mean one in which the main diagonal 
is comprised of sub-blocks for which each incremental sub-block contains one less row and column compared to 
the previous. In constructing these models we consider hierarchical Bayesian models with hyperparameters on the 
priors for development factors and specially developed matrix-variate priors on the covariance structures which 



161 and Merz and 



preserves the conjugacy properties of the independence models developed in Merz and Wiithrich 
Wuthrich [li- 

For the independent and Gaussian copula based PIC models we develop a class of conjugate posterior models 
that can be efficiently estimated via an MCMC sampler known as a block Gibbs sampler. However, the extension 
to general copula dependence structures requires non-conjugate priors, making it necessary to develop adaptive 
MCMC algorithms. Adaptive sampling uses previous iterates of the Markov chain to form more efficient Metropo- 
lis proposals for the parameters, this class of MCMC algorithm has received growing attention in the statistics 
literature since it was recently developed and is now recognized as an important tool for Bayesian inference. There 
is an increasing interest in utilizing adaptive MCMC to facilitate more efficient sampling ( Andrieu and Thoms , 
Atchade and Rosenthal [^). The adaptive techniques that we adopt in this paper fall within the general framework 



of adaptive Metropolis, and employ the optimal scale factors (Roberts and Rosenthal [25|) from the Single Com- 
ponent Adaptive Metropolis (SCAM) algorithm (Haario et al. [l^). There have been some initial utilisations of 
adaptive MCMC specifically in financial modelling such as [2l| and the references therein. In addition the adaption 
strategies we consider in this paper involve extensions of Euclidean space Adaptive Metropolis to the space of 
positive definite matrices, creating a class of matrix variate Markov chain adaptive proposals. 

In the mixture copula based PIC models, we design data augmentation strategies which are a class of auxiliary 
variable methods. We modify these approaches to the PIC copula based models in order to circumvent the challenge 
of intractable likelihood evaluations which arise form the structure of the PIC reserving triangle. In particular we 
argue that the tail dependence features of the model should be consistent accross all development years for both 
payment and incurred loss data. This poses an evaluation challenge for the likelihood as it involves evaluation of 
marginal likelihood quantities given the observed data in accident year i, given by payment and incurred losses. 
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The integral required when utilising mixture copula structures over the accident years is intractable, therefore we 
introduce auxiliary variables into the Baycsian model in a data-augmentation structure to overcome this dificulty. 

2. Review of the Merz-Wuethrich Independence Copula Paid-Incurred Claims Model 

This section introduces the PIC model which involves two sources of information. The first is the claims payment 
data, which involves payments made for reported claims. The second source of data incorporated into the statistical 
estimation are the incurred losses corresponding to the reported claim amounts. The differences between the 
incurred losses and the claim payments are known as the case estimates for the reported claims which should be 
equal once a claim is settled. This imposes a set of constraints on any statistical model developed to incorporate 
each of these sources of data into the parameter estimation. We use the constraints proposed in Mcrz and Wiithrich 
[igI ] which are used to specify a model based on a claims triangle constructed from vertical columns corresponding 
to development years of claims and rows corresponding to accident years. This structure for the observed data is 
summarized in triangular form which is utlised for both the claims payments and the incurred losses, as presented 
in Figure [TJ 

Without loss of generality, we assume an equivalent number J of accident years and development years. Further- 
more, we assume that all claims are settled after the J-th development year. Let Pij be the cumulative claims 
payments in accident year i after j development periods and lij the corresponding incurred losses. Moreover, for 
the ultimate loss we assume the constraint discussed on the case estimates corresponds to the observation that 
predicted claims should satisfy Pi j — li^j with probability 1, which means that ultimately (at time J) the claims 
reach the same value and therefore satisfy the required constraint. 

We define (i) Po-.j^-.j = {Pk.i : < k < J,0 < I < j}. (ii) Let A and B be square matrices. Then diag{A, B) is the 
diagonal matrix, with the diagonal elements of A appearing topmost, then the diagonal elements of B. Let the 
matrices A and B be as in (ii). Then the direct sum of A and B, written as A® B is the block diagonal matrix 
with A in the top left corner and B in the bottom right corner. It is clear that the definitions in (ii) and (iii) can 
be iterated. That is diag{A, B, C) = diag{diag{A, B),C) and A® B <3) C ^ {A® B) Q C. (iv) Define the d x d 
diagonal square identity matrix according to 1^. (v) Define the indicator of an event by the dirac-delta function 5i. 
(vi) Define the vectorization operator on a p x n matrix A, denoted by Vec{A) , as the stacking of the columns to 
create a vector. 

As in Merz and Wiithrich [l^ j we consider a Log-Normal PIC model as this facilitates comparison between existing 
results and the results we derive based on different dependence frameworks in extensions to this model. 

We now introduce the PIC model and the statistical assumptions for the independent case, followed by remarks 
on the resulting marginal posterior models for the payment and incurred losses. 

Model Assumptions 2.1 (Independent PIC Log-Normal {Model I)). The model assumptions for the in- 
dependent model of Merz and Wiithrich fldl l are: 

• The cumulative payments Pij are given by the forward recursion 

Pi ■ 

Pj.o = exp(C,.o) and =cxp(^,,j) forj = l,...,J 
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Figure 1: Claims triangle forpayment data and incurred data, ineluding constraint on zero case estimates at development period J 
(source Merz and Wiithrich llal V 



The incurred losses lij are given by the backward recursion 

h,j = Pi, J and ^^"^ = exp (-Cij-i) ■ 



The random vector (^cOi 



1 Cj,Ji Co,Oi ■ • ■ 1 Cj../-i) It-o-s independent components with 



for i G {0, . . . , J} and j G {0, . . . , J}, 
for k e {0, . . . , J} and Z G {0, . . . , J - 1}; 



• The parameter vector for the model is & ^ ($0j • • ■ j ^'o: • ■ • : ^J-ii cto, • • • , fj, tq, . . . , tj_i) . It is assumed 
that the components of are independent apriori. The prior density for © has independent components, 
with CTj, Tj both positive for all j . 



• It follows that 



log ( ^ N , a]) and log ( 



(2.1) 



Let {P, 1} = {Pi,j, Ik,i', l£ 'i,j,k,< J, < Z < J — 1}. Then, based on Model Assumptions 2.1 and the observed 
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matrices P and /, the likelihood for is given by three components, see derivation in Merz and Wiithrich 
Section 3.3, Equation 3.5]. The first and third components correspond to the payment and incurred data and the 
second component corresponds to the imposition of the restriction that ultimate claims for payments Pi j match 
incurred li^j for all accident years, giving: 



/(p,ii0) = nn exp|--i^(ci>,-iog(-^))^ 



n 



Component 1: payment 

1 



— CXp < — — ^ 5 



:{i^j-t - - iog( p'"^ ' )? 



Component2: Discounted final development year restricted payment and incurred 



Components: incurred 



where V^. = T.ir=0 C^m + H'LJ, ^? = Em=0 m = Em=0 ^ud [Ij Y:L=^ - ^IJ' 



(2.2) 
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As noted in Merz and Wiithrich [3] : there are several consequences of the model assumptions made regarding the 
restriction 7^. j = P^. j which applies for all i g {1, 2 . . . , J}. The first is that this condition is sufficient to guarantee 
that the ultimate loss will coincide for both claims payments and incurred loss data. The second is that this model 
assumes that there is no tail development factor beyond the ultimate year J. However this could be incorporated 
into such models, see Merz and Wiithrich [2]. 



Merz and Wiithrich discuss the relationship between the proposed Independent Log-Normal PIC model and 
existing models in the literature for payment loss based reserving and incurred loss based reserving. In particular, 
Merz and Wiithrich [2 [Section 2.1 and 2.2] show that the resulting cumulative payments Pi,j, conditional on 
model parameters 0, will satisfy the model proposed in Hertig [12[ and the incurred losses lij, conditional on 
model parameters 0, will satisfy the model proposed in Gogol 8|. Lemma 12.21 summarizes their results. 



Lemma 2.2. The relationships between consecutive payment development year losses in a given accident year is 
given conditionally according to 



log 



P 
P 



) P0:Jfl:j-l,Q - M a|) , Vj > (2.3) 

in agreement with Hertig 's model. With conditional moments given according to the Chain Ladder property as in 
Merz and Wiithrich fldl . Lemma 5.2] by, 

E [P,jPn:.LO:j-i, 0] = P^J~l cxp ($, + /2) . (2.4) 

Furthermore, conditional upon the model parameters &, for allO<j<j + l<J the relationships between 
consecutive incurred losses in a given accident year are given in Merz and Wiithrich fldl l [Proposition 2.2] according 
to 

[log (/,,,+0| /o:XO:,-i, 0] ^ + ^ (log(/..j) - > ^1+1 (1 - '^l+i/'^D^ , (2-5) 

These results are consistent with the model assumptions of Gogol, and are derived using properties of multivariate 
normal distribution, see Lemma 2.1 in Merz and Wiithrich tla] . 

Furthermore, for all accident years i € {1,2,..., J}, the resulting conditional expected ultimate payment loss 
equals the expected ultimate incurred loss, given the model parameters 0, and is expanded in terms of the model 
parameters according to Equation (|2.6p . which are given by Merz and Wiithrich [l6|. Equation 1.1] as, 

(2.6) 



E [P,,j| @]^E [/,,j| 0] = exp [ ^ + a^/2 j . 

\m=0 / 



3. Incorporating the Gaussian Copula into Paid-Incurred-Claims Models 

This section discusses an important aspect of extending the original Log-Normal PIC model of Merz and Wiithrich 
[l^. In particular, when this model was developed in the independent setting it was observed by the authors that 
the assumption of conditional independence between ^^j- and (k,i for all i, j, k , I e {1, 2, . . . , J} was not necessarily 

241 ] discovered evidence for strong linear 



consistent with observations. In particular, they note that Quarg and Mack 
correlation between incurred and paid ratios. In Section 3.1 we explore in detail a different approach to incorporate 
dependence structures into the Log-Normal PIC model. Some aspects of the new approach have subsequently been 
proposed in the literature, while others are novel developments proposed in our article. We note that in developing 
the extended models, the convenient properties of conjugacy in the Bayesian framework, which aids estimation, is 
often lost. Hence, after presenting the models we develop efficient state of the art statistical estimation strategies 
based on adaptive MCMC. 
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3.1. Dependence via Payment Loss Ratios and Incurred Loss Ratios (Model II) 



This secti on K cncralizcs the model by Happ and Wuthrich 



Wuthrich 



11 



ich ill, 



which has a static covariance structure, see Happ and 



Figure 1.1]. We use a Bayesian approach, based on results in Lemma FA . 2 1 and Model Assumptions 13 . 1 [ 
to estimate the extended models. We use properties of the matrix- variate Wishart and Inverse Wishart distributions 
to develop a Gaussian copula based statistical model. The relevant matrix-variate distributional assumptions and 
properties are provided in Lemma IA.2I and Lemma IA.3I 

Model Assumptions 3.1 (Dependent Payment-Incured Ratios: PIC Log-Normal (Model II)). The model 
assumptions for the Gaussian copula PIC Log-Normal model involve: 

• The random matrix G K(2J+i)x(2,/+i) j-Qpf^g^jiUjig ^/jg covariance structure for the random vector con- 
structed from log payment ratios {^i.j ~ log ("F^^^)) '^^'^ ^^9 ^i^curred loss ratios (^Ci,j = log ( j^' +i )) 

i-th development year, denoted by 'E^i = (^i, Oj Ci,i: Ci. 2, Ci, 2, Ci.j); assumed distributed according 

to an inverse Wishart distribution prior (see definition and properties in Lemma \A.2\ and Lemma [A.3\) . 

-ZW(A„fcO (3.1) 

where Ai is a ((2 J + 1) x (2J + 1)) positive definite matrix and ki > 2J. 

• Conditionally, given & = ($o, • • ■ , ^J, ^Oi ■ • • i ^j) the (2 J+ 1) x (2 J+ I) -dimensional covariance matrix 
E, we have: 

— The random matrix, constructed from log payment ratios = log ("F^^^)) ^'^'^ W incurred loss ratios 
,jj = \og(j^i^^\ denoted by E and comprised of columns Ei = (C,,o, 0,1, Ci,2, Ci,2, • ■ • , ■Ci.j, Ci.j), 



(c. 



is assumed distributed according to a matrix-variate Gaussian distribution /^^^^ (E!|M, E, see the 
definition and properties in Lemma \A.l\ The sufficient matrices are then the ((2 J + 1) x ( J + 1)) mean 
matrix M = [&',...,&'], column dependence given by ((2J + 1) x (2 J + 1)) dimensional covariance 
matrix E and row dependence given by (( J + 1) x ( J + 1)) dimensional matrix Q. If Q = Ij+i, the 
covariance of the vectorization of E = Vec{E) is 



^ J 

i=0 



/ Eo ... 
' El ... 



(3.2) 



V 







where it is assumed in the model in Happ and Wuthrich \1m J that E^ = Cow(E!i) = E. However, this 
need not be the case and it is possible to consider two extensions, the first in which Coi'(Si) varied as 
a function of i £ {0, 1, . . . , J} and the second being the most general of these model structures given by 
the assumption 



'(§ 



(3.3) 



— For all accident years, i G {0, 1, . . . , J}, the ultimate payment losses and incurred losses are equal a.s., 

• The matrix E is positive definite and the components of © are independent with prior distributions 

- (0„ si) and ^ M {,}j,,t]) , (3.4) 

and hyper-prior distributions 

sf ^ IQ {ai, (ii) and i| ~ IQ {aj,bj) , (3.5) 
for all i €z {1, . . . , J} and j G {0, . . . , J}. 

This model extends the model developed in Happ and Wuthrich \v\ which assumes that E is fixed and known with 
a tri-diagonal structure. The extension we consider generalizes the dependence structure to be unknown a priori 
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and given an inverse Wishart prior for matrix S, so it forms part of the inference given the data, in the Bayesian 

where they assume E = E;, Vi G {0, 1, . . . , J}, wo also 



inference. In addition, unUkc in Happ and Wuthrich 
allow for variation in E^ across development years. 



Given these model assumptions, we now consider two consequences of the proposed model structures for the 
dependence between the log payment ratios and the log incurred loss ratios given in Lemma 13.21 and Lemma 13.51 

Lemma 3.2. Conditional upon and ki, for all i in {0, 1, . . . , J}, and given the marginal distributions for E^ 
follow Ei ~ IW (Ai, ki) with Ai a ((2 J + 1) x (2J + 1)) positive definite matrix and ki > 2J, the joint distribution 
for the ((2J'^ + 3J+1) x (2J^ + 3J+1)) covariance matrix E for the vectorized matrix for S, given by 'E. = 
T^ec(S), under the assumption of independendence between development years, 

J 

S) =0E, = (Eo®---®Ej), (3.6) 



results in a joint distribution given by: 

E-Iw(A,fc), (3.7) 
with degrees of freedom k = X]/=o '^•^^ + ^"^ '^^'^ scale matrix 



A = 0A,. (3.8) 



i=0 

Furthermore, the joint prior mean and mode for the distribution of the random matrix A are 

1 



E 



E|A,fc 



A, and 



m ( E ) = ^- -. A. 

2 J2 + 3 J + 1 + Y.L^ 



(3.9) 



The proof of this result is a consequence of the results in Lemma \A.2[ the model assumptions and the properties 
of an inverse Wishart distributions; see Gupta and Nagar \^ [Chapter 3]. 

□ 

Remarks 3.3. We can demonstrate that under the proposed model assumptions the selection of the factorized 
covariance structure in Lemma \ 3.2\ produces Bayesian conjugacy in the joint posterior of the model parameters 
given observed payment losses and incurred losses. 

Remarks 3.4. It is noted in Happ and Wuthrich fTil / and Lemma \3.2\ that in formulating the likelihood structure 
for this dependent model it is more convenient to work with the one-to-one (invertible) transformation for the 
observed data defined marginally for the i-th development year according to 

[X,\&] = [BiEi\@] ~ M {B,AU,B,^,B[) , (3.10) 

where Mi is the i-th column of matrix M and Xi £ 1^2'^+-^ defined by 

Xi = [log /i,o, log P,;. 0, log /i^i, log log log log Ji^j] . This results in the joint matrix variate 

Normal distribution for random matrix X = [Xg, X'l, . . . , Xj] of all observed losses for payment and incurred data 
given after vectorisation X ~ Vec {X) by 



x\e 



BEl® 



■N{BVec{M),B{Y.®n)B'^) . (3.11) 



Furthermore, if we consider the property of multivariate Gaussian distributions given in Lemma 13.51 we can find 
for the i-th accident year the required conditional distribution of the unobserved claims for payment and incurred 
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loss data under the specified model. Furthermore, we can find the conditional distribution for unobserved claims 
for payment and incurred losses in the i-th accident year, given all observed claims triangles for payments and 
incurred losses data, see Lemma [3.51 below. This is directly relevant for specifying the resulting likelihood model. 



Lemma 3.5. Consider a (n x 1) random vector Y with multivariate Gaussian distribution, Y ~ A/'(/x, S) , where 



fi = [fii, . . . , fJ-n] o,nd Cov{Y) = E, and partition Y = 
given Y^^^ and the marginal distribution of Y^^"^ is 

y(i)|y^(2) 



^y(2) _ Then the conditional distribution of Y^^"^ 
-AA(/i,S), (3.12) 



with fi ~ fii + El. 2 {Y^"^^ — /^'■^^) ci'n'd the Schur complement E = Ei.i — Si,2S2 2^2,1 under the partitioning of 
the mean and covariance given by 



Ml 


and ^ ~ ^ 







El, 2 E29 



(3.13) 



Definition 13.61 below defines a family of permutation matrix operators. This permutation family allows the rep- 
resentation of the vectorization of the two loss triangles under different permutations that facilitate dependence 
specifications in the proposed models that admit conjugacy. 

Definition 3.6. Let Y be an n x n matrix, with Y = [Yi^i, ^1.2, • ■ •]' and with Vec{Y) defined as 
Vec{Y) — [Yi^i, Y1.2, . • . , Yi^n, ^2.1, • ■ • 7 ^2,71, ■ • ■ , Yn^n]' ■ Define the family of permutation matrix operators, denoted 
by and indexed by p x 2, p < , indices matrix (vector of tuple elements) i with j-th element [i ]j = 
{(fc, l);k,l & {1, 2, . . . , ri}}, and defined according to the mapping : Vec{Y) 1— >■ Vec{Y)* given by 



VI {Vec{Y)) = Vec{Y) 



(3.14) 



where we define Yjj j as the element of matrix Y corresponding to the resulting tuple index location in the j-th 
element (column) of (tuple vector) i , P* an x permutation matrix defined by 



Pi =Pi ei„2_ 



Pi 



0. 



71'^ —PjU-^ —p 



•^n-^ —p,n^ —p 
I 2 

■^71^ ~p 



(3.15) 



and Pi is a matrix with only non-zero identity elements at the p locations in the indices matrix tuples in i 
corresponding to index elements. 



Using the property of the multivariate Gaussian distribution in Lemma 1211 one can state the result in Proposition 
3.7l which is based on a generalization of the result in Happ and Wuthrich 11[ [Lemma 2.1] to the model developed 
above. We consider two cases for the dependence structures in Proposition 13 . 71 and Proposition [ 



Proposition 3.7. Consider the i-th accident year. Conditional on the model parameters and the covariance 
matrix of the i-th accident year 

[^i]l,l [^■i]24 
[^']l,2 [^i]2.2 



(3.16) 



the dependence structure 51 = Ij+i and the observed payment losses and incurred losses in the i-th accident year, 
denoted by X^^-* = [log/i,o,logP,;.o,log/i,i,logPi_i, ... ,log/i, log with Xi e M'', the conditional distri- 
bution for the log of the unobserved payment losses and incurred losses 
(2) 

{XI = [log/i,j_i+i,logPi,j_i+i,...,log/i,j_i,logPi,j_i,log/i,j]) is given by 



I. -(2) (2) 

where fi} — fi] 



XW\& 



[E,]2,,[E,]-;(xf)-/.f)) andE"/') 



A/-(/i(^\E-/^' 

= [^i]22- 



(3.17) 
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Proposition 3.8 (Conditional Distribution of Unobserved Payment and Incurred Losses). Consider the 
i-th accident year and define indices for this year (vector of tuples), given by matrix 

i = {(fc, j) : Vj e {J — A; + 1, . . . , J}} U {{k, j) : \fk G {0, 1, . . . , J} , j € {0, . . . , J — k}}. Then consider the trans- 
formed vector of log payment and incurred losses ( X ) defined by 



(3.18) 

correspond to all 



VI [X]^N [P* Vec{M),P: (E ® n) (Pi )') , 

(1) 



VI [X 



l:.7-i-l 



for which the first J — i elements of the permuted random vector 
un-observed payment and incurred loss random variables, and the remaining J — i to J — i + (^2n=-ii:^ ^ '^)) 



(2) 



Th 



en, 



elements are the observed payment and incurred data, denoted X* = Vt \ X , 

^ ./-*:./-*+(E„'=-i(./-")) 

conditional on the model parameters 0, the general dependence structre S = E 17 with matrices S and 17, and 



X' 



(2) 

the following results hold: 

The conditional distribution for the log of the unobserved payment losses and incurred losses in the i-th year, 

X* = I A I IS given 



corresponding to the first J — i elements of the permuted random vector 
by 



vt (x) 



X* 



(1) 



X* 



(2) 



,0 



(1) 



(3.19) 



• The covariance matrix e/"'^^ is the postive definite 

- n))) X (j - z + (E,t-i(^ - " 
top sublock of the permuted covariance matrix 

p* (E <E> n) {p* y = 



sub-matrix denoted below by T and defined by the 



r [p: {E^Q)iPn'],, 

[p* (s ® n) (p* )'] ^ 2 [p: (s ® Q) [Pi y] ^ ^ 



(3.20) 



• Given, this covariance matrix one specifies the conditional mean vector, denoted by 

(2) 



X' 



spect to the first J—i elements 



X* 



(1) 



/i*-^^ I , according to the subblocks of the T covariance matrix defined with re- 



and remaining elements of 



X* 



(2) 



dl asu^^^ = [p.* Vec{M)]^^ 



J-i 



and the second J^i to J~i-\-(^J2n=-ii'^ ~ dements are given by fi^^^ = [P^* yec(A/)] ^ i+{J2'' iJ "))' 

Having specified these statistical assumptions, we can formulate the joint likelihood from the observed data for both 
payments and incurred claims conditional upon the model parameters according to Equation p.2ip . The joint data 
likelihood function in the dependent Log-Normal PIC Model I for the random vector of observations corresponding 
to the first E^^^_j^( J — n) elements of the permuted random vector, given by X* = VI \X] , 

where we define indices in this case by i = {(«, j) : Vi G {0, 1, . . . , J} , j G {0, . . . , J — i}}. The resulting likelihood 
is given by the matrix- variate Gaussian distribution in Equation p.2ip . 

1(1) 



exp 



©,S,f7 1 = 

(1) 



[X*]*" - [V* {Vec{M))\''''>\ [[P* (E®n)(P;)f ' {\X*Y^ - [vt (V^ec(M))]'^ 



(3.21) 



We note that our proposed models also allow one to consider the dependence structures of Happ and Wuthrich 
[ll| who assume that — S, Vi e {0, 1, ... , J} and Vt ~ Ij+i, with the specific setting of E via a tri-diagonal 
correlation matrix with three correlation parameters which are assumed either known a priori or estimated prior to 
inference in the PIC model. Such an approach was motivated by the belief that a positive change in incurred loss 
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results in an immediate payment in the same development period, and the remaining increased expectation is paid 
with some settlement delay. Therefore, the incurred losses increments C| arc assumed to be positively correlated to 
the claims payments increments ^ij, Cij+i ^i,j+2 with positive correlations po,pi,P2, respectively. However, 
the argument for more general dependence structure that were introduced as extensions to the model of Happ and 
Wuthrich [ll| are developed to account for the fact that these assumption may not be true, especially in long tail 
portfolios, such as compulsory third party. If the status of a claimant changes and requires long term medical 
treatment and rehabilitation, it might result in substantially high loss in the subsequent lengthy lag periods. The 
paper also assumes that the dependence is the same across different lag years, which is not always a realistic 
feature of such data. Our article aims to fill this gap and enhance the correlation structure in PIC models whilst 
maintaining a parsimonious model specification. 



3.2. Dependence Between Development Lag Years for Payment Losses and Lncurred Losses (Model III) 

This section considers an alternative dependence structure motivated by the fact that dependence between lag years 
is practically appealing in claims reserving practice. It affects the estimation of outstanding claims the most, and is 
widely recognized by actuaries in claims reserving. Lag is the measure of the difference between incurred month and 
paid month. Depending on the nature of the portfolio, many insurance claims often have lengthy settlement periods 
due to various reasons such as late reported claims, judicial proceedings, or schedules of benefits for employer's 
liability claims. During the lengthy lag periods, large payments in the previous lag period normally follow by small 
payments in the subsequent lag period. There may in fact be positive correlation if all periods are equally impacted 
by a change in claims status, e.g. if a claim becomes litigated, resulting in a huge increase in claims cost. There 
may also be negative correlation if a large settlement in one period replaces a stream of payments in later periods. 
The model developed in this section mainly focuses on capturing this feature of dependence between lag years. To 
achieve this we propose a block covariance structure for the covariance matrix, which will respect the dependence 
between each lag point. The model we propose is summarised in Model Assumptions 13.91 below. 

Model Assumptions 3.9 (Dependent Development Lag Years: PIC Log-Normal {Model III))- The fol- 
lowing statistical model assumptions are developed: 

• Let Sf e §D+(J - i) be the {J - i) x {J ~ i) random covariance matrix on the space SD+(J — i) of positive 
definite covariance matrices of dimension {J i) x (J — i) corresponding to the observed payment data 
[log Pi fi, log Pi ^i, . . . ,logi-i^j_,] in the i-th accident year and analogously for incurred loss data G §B+(J — 
i). When i = we consider G §B+(J + 1) and for incurred loss data log Jo,o:J-i with Sq G §D+(J). 
Assume an inverse Wishart distribution (see Lemma ] A. 3\ and Lemma \A.S\) for each matrix defined according 
to 

Sf ^ IW (Af , fcf ) and Sf ~ IW (Af , fcf) , (3.22) 

where Af and A^ are the inverse scale matrices for the prior for the payment and incurred loss data covariance 
priors respectively. Hence, the joint covariance between all observed payment and incurred loss data satisfies 
the telescoping diagonal block size covariance structure: 

E = C ( [log Po,0 , ■ • ■ , log Po,J , log Pl^O , log Pl,J- 1 , . . . , log P,7,o , log /o,o , ■ • • , log /o, J- 1 , • • ■ , log Ijfi] ) 

= f e So") © f e ^i] - f f e ao-) e f a^) , x: (fc.r + h) 

\i=0 / \i=0 / \\j=0 / \i=0 / i=0 

• Conditionally, given = {^o, . . . , $j, ^o, . ■ . , ^j) and the covariance matrix S, we have the following results 
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— Consider the marginal distribution of the first ^X!n=-i("^ ^ elements of the vectorized random matrix 
of observed payment and incurred losses, with i-th column Xi G given by 

Xi = [log /i,o , log Pi,0 , log /i, 1 , log Pj, 1 , . . . , log j_ 1 , log Pj, j_ 1 , log /i, j] . 

Then given the matrix of permutation indices i ~ [(1, 2), (1, 4), ... , (1, 2( J — 1)) 

,(2,2), (2,4),..., (2, 2J- 4),..., (J, 1), (1,1), (1,3),... (J - 1,1), (J - 1,2)] characterizing the elements 
of the marginal distribution for the observations, the transformV^ {Vec{X)) has multivariate Gaussian 
distribution with covariance structure E. Note, {Vec{X)) = [log Pq.o? log Po.i: ■ • ■ jlogPo.J 
, ■ • ■ ,logPxo,log/o,o,- ■ • ,log/o,j-i,log/i,o, • • ■ ,log/,/-i,o,log/j-i,i]. 

— For all accident years, i g {0, 1, . . . , J}, the ultimate payment losses and incurred losses are equal almost 
surely, Pij = h^j. 

• The matrix S is positive definite and the components of are independent with prior distributions 

$^ - A/" (0„ sj) and *j - N [il^j^t]) (3.24) 

and hyper-prior distributions 

sf ^ ig (a„ A) and t] ^ TQ {a,,b,) (3.25) 
for all i & {1, . . . , J} and j e {0, . . . , J}. 

This proposed model is therefore another generalization of the dependence structure of the model structure proposed 
in Happ and Wuthrich [ll| . As such, the likelihood structure is given by the multivariate Gaussian given in Equation 
(|33T|) with the covariance matrix given by the telescoping diagonal block size covariance matrix structure in 
Equation ([X^ . 



3.3. Hierarchical Bayesian Conjugacy Under Gausian Copula Dependent PIC: Models I, II, III 

Under the Gaussian copula based dependence models, the ability to obtain the observed data likelihood in the form 
of a multivariate Gaussian distribution means that we obtain conjugacy properties. This makes the estimation of 
such models by MCMC more efficient because we can us Gibbs sampling in blocks. This section presents a generic 
set of such conjugate models for any of the dependence structures specified in Models I, II and III. 

Lemma 3.10. Conditional upon the parameters & and the covariance matrix E, the permuted data VI {Vec{X)) 
can be transformed to produce the independent likelihood in Equation ^2.2^) . This is achieved by considering the class 
of vector transformations T : R'-''^^-' i— >■ M'^''^^-', such that if the initial covariance structure of random vector X 
was given byTi = Cov{X), then the resulting covariance structure Cov{T{X)) = Id. The required rotation- dilation 
transformation is obtained by the spectral decomposition of the covariance according to a spectral decomposition (see 
Stoica and Moses ^2^1) S = UAiu' where U is a (d x d) matrix of eigenvectors of S and A is a diagonal d x d 
matrix of the eigenvalues ofT,. Therefore the following holds for each of the models under a transform of the vector 
of permuted observations T {V^ (V^ec(X))); 

1. Model II - When S = S ® 17, with fl = Ij+i then, T {V^ {Vec{X))) = (uA^ ® Ij+i) -p; {Vec{X)) , where 
the ((2J+ 1) X (2J+ 1)) covariance E is decomposed as UA^U'. 

2. Model II - When E 0/^^ S„ T (P? {Vec{X))) = U^Af^ Vt {Vec{X)) , where each accident year's 
dependence between payments and incurred losses is given by the (2J + 1) x (2.1 + 1) matrix which is 
decomposed as UiA^U[. 

3. Model III - When t = (©t^ S^) (©t^ S^) , 

r (P; iVeciX))) = (^0C/f (Af)^^ ® (^0L// {Aiy^ VI {Vec{X)) 
where each of the covariance matrices Sf and decomposed to Uf (Af ) ^ [Uf ]' and (Af) ^ {U^Y ■ 
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In each case, the resulting transformed random vector T (V^ {Vec{X))'^ , with elements Pij and Iij, will produce a 
likelihood model given for the transformed data according to the independent Model I of Merz and Wiithrich fldl l as 
defined in Equation h2.2]) . Of course this is defined now with respect to components in the likelihood corresponding 
to the transformed components, as detailed in Equation \3.11\) . 



Remarks 3.11. The consequence is that results in Lemma \S. 10\ are that the conjugacy properties derived for the 
independent model in Merz and Wiithrich fldl l can be directly applied post-transformation. This is of direct interest 
for MCMC based sampling schemes. 

In the models described so far, the fohowing full conditional posterior distributions are now of relevance to the 
Bayesian MCMC estimation procedures developed for Models I, II and III. 



Lemma 3.12. The full conditional posterior distributions for sub-blocks of the model parameters can be decomposed 
under Model I, II and III into a conjugate model. 



• Conjugate Posterior Distribution for Development Factors: under the transformations T {Vl (X^ec(X))) 
on the data, described in Lemma \3.10\ the full conditional posterior distributions for sub-blocks of the trans- 
are given by (see Merz and Wiithrich flA l [Theorem 3.4] for the 



formed model parameters ^$o:,/i ^O:, 
independent case): 

$0:,/,*0:,/|S,f^,r(P; {Vec{X)))] ^Af{n,A) 

with posterior mean 11 and posterior covariance A, where the components of A^^ = (a,i,m)Q>„ ,„<2j 
given by 



(3.26) 

are each 



(n-l)A(?ri-l) 

ln,m = {Sn^ + {J - n + 1)ct^^) 6n=m + ^ 



'^f - ^i) \ for < n,m < J, 



aj+i+n,j+i+m = {tn + {J - n)T,^ ^) 5n=,n + ^ [vf - ) \ for < u, m < J - 1, 



(3.27) 



(n— l)Am 



0'n,J+l+r, 



A 



n.J-\-l-\-m 



{iyf~ujf) , /or < n < J,0 < m < J- 1; 



4=0 



where 5n=m is the indicator of the event that index m matches n, m An is the minimum of m and n and the 
posterior mean is given on the transformed scale by. 



*0:J,*0:,/ 



with 



Pr 



i=0 
.1-3-1 



. Pi J -I 



1=0 \^«J + l/ i=J-i 



A(co,ci,...,cj,bo,...,bj 



J 



(3.28) 



log 



(3.29) 



Given the transform vector 3>0:J!^0:J , the parameters on the orginal scale can be expressed according to 
the unique solution to the system of linear equations: 

1. Model II - On the untransformed scale, the solution is given by the following system of equations 

[$0:J,*0:,/]' = C/"'A-^ 



$0:J,*0:J 



(3.30) 



2. Model II - On the untransformed scale, the solution is given by the following system of equations for 
each i G {0, 1, . . . , J}, where we can randomly select i or deterministically scan through i for the results, 



[<i>0:J,*0:j]' = C/r'A; 



0:,/ 



(3.31) 
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3. Model in - On the untransformed scale, the solution is given by the following system of equations, 

J 



[$0:J, $0:J-1, *0:J-2, ■ • ■ , = {Uf) ' (A 

i=0 
J 



i=0 



$0:,/,$0:J-l,$0:J-2, 
*0:J,*0:J-l,^0:J-2, 



• Conjugate Posterior Distribution for the Covariance Matrix: Given the transformed observed pay- 
ment and incurred losses have a multivariate Gaussian likelihood, as specified in Equation h3.21\) . with co- 
varaince matrix S = E Cg) O and mean vector Vec{M). Then the posterior for the covariance matrix is the 
Inverse- Wishart- Gaussian distribution detailed in Peters and Mellen Fil / [Section 3] and Peters and Godsill 



S|$o:J,«'o:J,r(n* {Vec{X)))\ ^IW[K + T{Vl {V ec[X))) T {Vl [Vec{X))) ,dtm{Vec{X)) + k 

In cases in which the covariance matrix T, takes any of the block diagonal forms presented in Models II and 
III, we may utilise Lemma \A.2\ and the result in Equation \3.12\l to further decompose the posterior covariance 
into blockwise components. 

Conjugate Posterior Distribution for the Hyper-Parameters on Development Factors: For all i 

we have the following Inverse Gamma- Gaussian conjugacy for the hyper parameters in Models II and III, 



and [til'^i] 



ig iai + 



We next present alternative tail dependence structures for the PIC model. Previous studies on claims reserving 
that have incorporated copula based models, such as Zhang and Dukic [2^ have done so through regression based 
frameworks. Zhang and Dukic [2 develop a parametric copula model to account for dependence between various 
lines of insurance claims. Their paper considers a bivariate Gaussian copula model with marginal generalized 
linear models to capture the positive correlation between the two insurance lines. Our article significantly extends 
the dependence modelling capability of the PIC model structure remaining in the frameworks presented above. 
However, to do so requires the introduction of auxiliary variables to enable computation. The approach developed 
involves modifying the posterior distribution by cmbcding the target posterior distribution for the model parameters 
into a much higher dimensional support comprised of the original model parameters and the additional auxiliary 
variables. The reason for this expansion of the posterior dimensions will be come clear below and is in general 
known in Bayesian statistics as an auxiliary variable framework. 



4. Incorporating Mixture- Archimedean Copula Dependence Structures into Paid-Incurred-Claims 
Models: Model IV 

This section presents an alternative parameteric approach to modelling and capturing dependence and tail de- 
pendence in the PIC model structure which involves considering copula based models within the PIC reserving 
framework. The dependence can be considered over the following combinations such as: 

1. Independent accident years and dependence between payment losses over the development years; 

2. Independent accident years and dependence between incurred losses over the development years; 

3. Independent accident years and dependence jointly between payment and incurred losses over the development 
years via a mixture copula, hierarchical copula (HAC) as in Kurowicka and Joe or a vine copula (d-vine, 
canononical vine) e.g. Aas et al. [ll; 
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4. Dependent accident years and independent development years for payment, incurred or both sets of losses. 



Our article concentrates on the mixture copula model which allows for combinations of upper and lower tail 
dependence of different strengths. We detail the class of auxiliary variable methods known in statistics as Data 
Augmentation and demonstrate how this class of models can be combined into our modelling framework to allow 
for consistent use of copula models in the PIC framework. There are many variations that can be explored in this 
approach. We give one such approach for Model IV, Assumptions 14.31 that is directly comparable to that used for 
Model II in Assumption 13.11 



We present fundamental properties of members of the Archimedean family of copula that we consider when con- 
structing mixture copula models in the PIC framework in the Appendix, see Lemma IB. II for the characteristics of 
the Archimedean family of copulas and Lemma IB . 2 1 for the required distribution and densities for three members 
of this family. In addition references Dcnuit et al. 0],Aas et al. [l| ,Embrcchts [?], Min and Czado [l^] and Patton 



provide more detail. 



In Lemma IB. II the property of associativity of Archimedean copula models is particularly useful in the PIC model 
framework as it allows us to obtain analytic expressions for the likelihood structure of the matrix- variate PIC model. 
This is particularly useful if one specifies the model as a hierarchical Archimedean Copula (HAC) construction. 

We consider the following popular members of the Archimedean family of copula models, due to their analytic 
tractability, their non-zero tail dependence properties and their parsimonious parameterizations. In addition, gen- 
erating random variates from these class of models is trivial given the generator for the member of the Archimedean 
family of interest. Lemma IB. 21 in the appendix presents the three Archimedean copulas for Clayton, Gumbel 
and Frank copulas that we consider and their properties. We use the following notation for copula densities we 
consider on [0,1]'', see Nelsen [l^ Section 4.3, Table 4.1] and Lemma lB.21 the Clayton copula density is de- 
noted by (ui, ...,Un; p'^) with £ [0,oo) the dependence parameter; the Gumbel copula density is denoted 
by c'~^{ui, ...,Un', p^) with p'~^ e [l,oo) the dependence parameter; and the Frank copula density is denoted by 
c^{ui, ...,u„;p^) with p^ G M/jO} the dependence parameter. 

In addition, we also note that the properties of these copulas of interest include that the Clayton copula does not 
have upper tail dependence, however its lower tail dependence can be expressed as Al = 2^^/'''^. The Gumbel 
copula does not have lower tail dependence, however its upper tail dependence of the Gumbel copula can be 
expressed as Xu = 2 — 2}lf . The Frank copula does not have upper or lower tail dependence. 

In this class of copula dependence models we consider the marginal distribution of each log payment or log incurred 
loss as distributed according to a Gaussian distribution and the joint distribution vector is modelled via a mixture 
copula comprised of the above three components from the Archimedean family. Such a copula construction will 
still produce a copula as shown in Lemma [4. II 

Lemma 4.1. Consider copula distributional members Ci (ui, U2, ■ ■ ■ ,Un) € A" , where defines the space of all 
possible n-variate distributional members of the Archimedean family of copula models, specified in Lemma \B.2[ 
Any finite mixture distribution constructed from such copula components that admit tractable density functions 
Ci (ui, U2, . . . , Un), denoted by c (ui, U2, . . . , w„) = X]"=i '^i'^i ^2, • ■ • , w„), such that X^IIli "^i ~ ^' 
density of a copula distribution. 

The proof of Lemma 14.11 is provided in Appendix O 
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4-.1. Understanding Bayesian Data Augmentation 



The modeling framework of Data Augmentation in the Bayesian framework is typicaUy invoked to deal with 
situations in which the likelihood evaluation is intractable to perform point-wise. This would make Bayesian 
inference in such a model also generally intractable. For example if one considers the generic likelihood p (j/iml^) 
with observation random vectors l^im, which can be evaluate point- wise as a function of parameter vector 9 with 
respect to a realization of the observation process yi-.n- 

In the setting wc encounter in the PIC models, we can generically consider the data random vector observation is 
partitioned into two vector sub-components Y = [l^*-^-*, Y^^'>'\ , of which only one component, say Y'-^'f, is actually 
observed. Then evaluation of the likelihood pointwise for 9 given a realization of y}^^ would require solving the 
integral in Equation 14.11 

p(yI^\0) = Jp{Yl}MYl'j)p{Yl^h0) dYl^. (4.1) 
Generally, this integral will not admit a closed form solution. Therefore, the Bayesian Data Augmentation approach 
involves extending the target posterior p (j)\Y^l^^ which is intractable due to the intractability of the likelihood 
to a new posterior model on a higher dimensional space, in which the target distribution is a marginal as given in 
Equation |42] 

p (d, yI^* \yI}^) = 1 (4.2) 

where V^.„ are auxiliary random vectors with prior distribution p I ^i.„ 'augmented' to the posterior pa- 

rameter space to allow tractability of the posterior inference. This will be explained in detail for the PIC copula 
models below. 



4-.2. Data Augmentation in the Bayesian PIC Copula Models 
Definition 14.21 gives some useful notation for the results that follow. 



Definition 4.2 (Auxiliary Data for Data Augmentation). Consider the defined loss data under the one-to- 
one (invertible) transformation for the observed data given by the joint matrix for all observations and auxiliary 
variables given by X = [Xq, X(, . . . , X'j] . In this framework, the i-th accident year is defined according to, Xi = 
[log/i.o, log Pi.o, log logPi.i, . . . , log log , log/i. j]. Consider the permutation of each vector of log 

payments and log incurred losses given by 

Xi = V- (Xi) = [logPi^Ojlog/'j,!, . . . ,logfi,j,log/i,o, log/j,i, . . . , log/i,j_i] . Now consider the further partition 
by the decomposition of observed log payment losses and unobserved log payment losses as well as these quantities 
for the incurred losses defined for the i-th accident year by. 



V 

~ ^i.obsT ^i.auxT ^i.obsT ^i,av 



x^ 

^O,i,o6s5 ■ • ■ 1 J—i,i,obs^ J — i+l,i,aux: ■ 



' J.i.aux 1 



x^ 



; J—i,i,obs^ J—i+l,i,aux: ■ ■ • : J—l,i.au 



. . . ,logPi„/_i,logPi,j_i+i, . . . ,logPi,j,log/i,o, • ■ • ,log/j„7-i,log/i,,7_i+i, . . . ,log Ji,j_i 



observed Payments 



unobserved Payments 



observed Incurred 



unobserved Incurred 



(4.3) 



Note, the introduction in this section 



Therefore the total data matrix of losses is given by X 
of the notation subscripts obs and aux allows us to make explicit the fact that the upper triangle of log payment 
losses and the upper triangle of log incurred losses are un-observed quantities for these random variables, while 
the lower triangular regions for such losses are observed. We denote these random variables as auxiliary variables 
(augmented) to the observed data random variables to create a complete data set of all losses. 
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By considering the unobserved data in the lower payment and incurred loss triangles as auxiliary variables to 
be jointly estimated along with the model parameters, we will demonstrate below that only under this approach 
is consistency ensured in the copula structure of the PIC model. However, we first make the following model 
assumptions about the statistical features of the PIC model. 

The following assumptions illustrate a choice of copula models for the mixture from the Archimedean family. 
However, there are many related specifications and frameworks that can be explored in this context, be we leave 
that to future research. 

Model Assumptions 4.3 (Data- Augmented Mixture Copula PIC {Model IV)). The model assumptions 
and specifications for the copula model we develop involve: 



with 



• Let the random matrix G R(2'^+^)^(2./+i) covariance for Xi = Xf^^^, Xf^^^^, X/^^^, X/^^ 
Xi G M^"^+"'^ for all i = 0, . . . , J . We assume that S is diagonal where 

E,,,^ig(a,;,A), V^e {0,..., J}, (4.4) 
where Ui and Pi are the known hyper-parameters for shape and scale. 

• MARGINAL DISTRIBUTION: given = ($o, • • • , ^'o, ■ • ■ , *j) and covariance matrices G 
]g(2j+i)x(2j+i) ^g Qgg^jYie the marginal distribution of the random matrix, of all log payments and log 
incurred losses X , comprised of columns Xi for the i-th accident year isjmatrix-variate Gaussian with density, 
defined as in Lemma \A.1[ with the (2 J + 1) x ( J + 1) mean matrix M = [©', . . . , ©'], column dependence 
given by (2J + 1) x (2J + 1) covariance matrix S and row dependence given by ( J + 1) x ( J + 1) matrix f2. 
Here we only consider the case of Q = Ij+i for the marginal independent case. 

• DATA A UGMENTED PIC MIXTURE COPULA LIKELIHOOD: Given X^„„^ , X^^a«:. - ■ • : ^j--i,aux' 
Kau.,Kau., X'j-i,au.^ ® ^ (*o, • ■ ■ , *J^*o, ■ ■ ■ , covartancc matrices S, e m(2J+i)x (2.7+1) 

p, the joint distribution of the random matrix (X) of all log permuted payment and incurred losses is assumed 
(in this example) to he independent between accident years. For the i-th column (corresponding to i-th acci- 
dent year), the joint distribution of all losses {Xi) is assumed to be hierarchical Archimedean Copula (HAG) 
mixture copula specified by distribution, 

X 



^P. ( ^ ( ^i.obs^ ^i.auxT ^l.obsT ^l.aux'^ ^ j 

(4-5) 

Cjp [f (X^„,„ Xf_ • [M]f„ S) ) Cl^ [f (x/,„,,, X/,_; [M]i^, s) ) , 



with supper script P and I denote the components for the log payments and log incurred losses in the i-th 
development year respectively and the density is given by 



f( X^ X^ ITA/fl ■ y 

J \-^i,obs''-^i,aux^-^i,obs^^i,aux\l-'-*^-'-\*iT'^^h'i t Hi 



_ _ _ _ 2,7+1 ^ (4.6) 

- [f (x,^,„ X,f,„,; [M]^„ S)) d'^^ (f (x,^„,„ X,^_,; [M]U, e)) [] M,, „ E,,,), 

where 

+ (1 - - ii;2)c^(c,s, (xf,_„,,;Mi^„Si,i) (x|,^,„,; M|„ Ej^j) ) , Se{P,I}, 

and such that ifi + W2 + (1 ^^t^i — ^2) = 1. This specifies a mixture of central, upper and lower tail dependence 
as denoted by the mixture of Archimedian copula models made up of Frank, Glayton and Gumbel members, 
such that for the source of data S , the copula parameters for each Archimedian family member is given by 
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p(G,S) ^ pI^'^"^ > 1 and pl^'^^ € M/{0}. Therefore the total conditional distribution corresponding to the 
likelihood model considered is given by, 

i=0 



exp I — 5^?" 



Copula Dependence in Data Augmented PIC Likelihood 

(4.7) 



(2j2+3J+l)/2 |f^|(2,/+l)/2 |^|(,/+l)/2 



Marginal Distribution in Data Augmented Likelihood PIC Model 

• Assume that the tail dependence features of the Data- Augmented copula PIC model are such that the 
dence structure is homogeneous accross accident years, = pf and p^ = pi for a/H G {0, 1, 2, . . . , J}. 

• Conditional on S, $ = [$o, ^"1, • ■ • , $j] and ^ = [^a, 5*1, . . . , ^ j] the hierarchical prior distribution on the 
auxiliary payment data for the i-th accident year is given by a normal distribution, centered on the development 
year mean, 

X[^au. - $^.,+2, . . . , <i>j] , S^) ' (4-8) 

The hierarchical prior distribution on the auxiliary incurred loss data for the i-th accident year is given by 

Xlau. - ^{[^J-^+l,^J-^+2, *,/] , ^'2) , (4.9) 

with S2 the lower portion of covariance E corresponding to the lower triangle matrix from ( J — i + l) through 
to J for all i € {0, 1,2, . . ., J}. 

• For all accident years, i G {0,1,..., J}, the ultimate payment losses and incurred losses are equal a.s., 
Pi,j = h.h P- a.s. 

• The matrix E is positive definite and components of are independent with prior distributions 

~ M ((/.„ s2) and ^ M {^Pj,t^) (4.10) 

and hyper-prior distributions 

sf - ig {a^,Pi) and t) - IQ [a,,b,) (4.11) 
for all i £ {\, . . . , J} and j G {0, . . . , J}. 

• The matrix S is distributed as E ^ ZW (A, k) and the copula parameters are distributed as p^'^ ^ IQ {oP , /3'~^) , 
p^^P ^ ig {a^, p^) and p^'^ ^ M (O, a^) 

Hence, we have made precise the auxilhary data scheme used in formulating the Data-Augmented-PIC model. In 
particular illustrating the importance of the role of the auxiliary data in evaluation of the model and estimation of 
the PIC claim development factors. Also we note we get indirectly via the data augmentation the distribution for 
the predicted payment and incurred Loss reserves. 

Remarks 4.4. The following remarks provide motivation for the Data- Augmentation and resulting incorporation 
of auxiliary payment and incurred Losses data. 

• The use of data augmentation in the above model structure is critical in the PIC model formulation, since 
it allows one to ensure that the dependence structure considered (in this case a HAC-Mixture) is consistent 
both across accident years and across development years. 

Note: In the case of a linear dependence structure such as with a covariance / correlation matrix under a 
Caussian Copula or Independent Copula model, such as those presented previously under Models I, II, III, we 
have that conditional distributions and marginal distributions are Caussian. This means that the evaluation 
of the likelihood is analytic without the need for auxiliary variables. 
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• In order to evaluate the likelihood one has two choices, to evaluate the observed data likelihood (Equation 
^.12^ ) or to evaluate the full data likelihood (Equation J^.Tp j. 

— The PIC copula model equivalent of Equation \4-2\ is the observed data likelihood is given for the i-th 
accident year by 

P [-^[obs, ^lobsl®, S> p) 

= j - jp {xlobs.xU\@, s, n, p, x^_) p x^_|0, s, 17, p) dxl^jxl^^ 

= /•••/ c^P (f [M]f„ e)) ~c'^, [f (x^,,„ X^,„,; [M]i„ s)) 

X -^i^^f^ {.^i,aux] Mi,aux, © S2) f^^^ (^ii -^'^ii ^) '^-^i^aiix ^Xj 

where matrix-variate Gaussian distributions fx^^ (^^d fx^^ o-fe as defined in Lemma \A.l\ 

with Xi^aux = Vec {Xf,^^^, Xl^^^ , Mi^aux = V^ec([$J_,;+i:j]', 

Xi = [^fofcs' ^z>u2;. ^i^ofcs' ^z>u:r] ^^c^ = [*o, ■ • • , $j , ^-Q , • • ■ , ^/le equivalent mean. 

— Clearly, the marginalization required to evaluate the Observed data likelihood involves intractable inte- 
gration, except in special cases in which the copula models are Caussian or independence copulas. 

• The full data likelihood comprised of observed and auxiliary data involves incorporating auxiliary variables to 
represent the unobserved data in the lower reserve triangle for payment and incurred loss triangles. These 
become part of the inference procedure and are required to be estimated jointly with the model parameters in 
the estimation methodology. 



5. Estimation via Adaptive Data- Augmented MCMC for Claims Reserving PIC Models 

It has been shown for the Independent and Gaussian copula models that we can obtain the observed data likelihood 
analytically. Therefore the posterior distribution for all the model parameters can be sampled via a MCMC proce- 
dure comprised of block Gibbs sampler updates. In the case of a more general copula dependence model in which 
the observed data likelihood cannot be analytically evaluated pointwise, we must resort to a Data Augmentation 
scheme. In this case we will be able to perform sampling via a general MCMC Metropolis-Hastings sampler. In 
particular we will consider automating such a sampler using an adaptive MCMC scheme. 



5.0.1. Adaptive Metropolis within Data- Augmented Copula PIC Models 

This section presents the adaptive proposal we use to sample the parameters and the auxiliary variables. The 
advantage of an adaptive MCMC mechanism is that it automates the proposal design through consideration of a 
proposal distribution that learns the regions in which the posterior distribution for the static model parameters 
and auxiliary data has most mass. As such, the probability of acceptance under such an on-line adaptive proposal 
is likely to improve as the iterations progress and the generated MCMC samples will ideally have reduced autocor- 
relation. In such cases the variance of Monte Carlo estimators of integrals of smooth functionals formed from such 
samples will be reduced. 

n 

There are several classes of adaptive MCMC algorithms, see Roberts and Rosenthal |25j. The distinguishing feature 
of adaptive MCMC algorithms, compared to standard MCMC, is the generation of the Markov chain via a sequence 
of transition kernels. Adaptive algorithms utilize a combination of time or state inhomogeneous proposal kernels. 
Each proposal in the sequence is allowed to depend on the past history of the Markov chain generated, resulting 
in many possible variants. 
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Haario et al. [10| develop an adaptive Metropolis algorithm with proposal covariance adapted to the history of the 
Markov chain was developed. Andricu and Thorns [31 is presenting a tutorial discussion of the proof of ergodicity of 
adaptive MCMC under simpler conditions known as Diminishing Adaptation and Bounded Convergence. We note 
that when using inhomogcneous Markov kernels it is particularly important to ensure that the generated Markov 
chain is crgodic, with the appropriate stationary distribution. Two conditions ensuring ergodicity of adaptive 
MCMC are known as Diminishing Adaptation and Bounded Convergence. These two conditions are summarised 
by the following two results for generic Adaptive MCMC strategics on a parameter vector 6. As in Roberts 
and Rosenthal [25[, we assume that each fixed MCMC kernel Q^, in the sequence of adaptions, has stationary 
distribution P {■) which corresponds to the marginal posterior of the static parameters. Define the convergence 
time for kernel Q.y when starting from a state 6 E E, as Af^ (9, 7) = inf{s > 1 : {6; — P (•) || < e. Under these 
assumptions, they give the following two conditions which are sufficient to guarantee that the sampler produces 
draws from the posterior distribution as the number of iterates tend to infinity. The two sufficient conditions are: 

• Diminishing Adaptation: lim„„).ooSupgg£;||(5rs+i (^j ') ^Qr^ (^j ') \\tv = in probability. Note, are random 
indices. 

• Bounded Convergence: For e > 0, the sequence {Mg {6,Tj)}°°^Q is bounded in probability. 
The sampler converges asymptotically in two senses, 

• Asymptotic convergence: limj_).oo||'Caw (9) — P (9) \\tv = in probability. 

• Weak Law of Large Numbers: limj_>ooi J2i=i 4>{(^) ~ J (l>{9)P{d9) for all bounded (j) : E ^ R. 

In general, it is non-trivial to develop adaption schemes which can be verified to satisfy these two conditions. In 
this paper we use the adaptive MCMC algorithm to learn the proposal distribution for the static parameters in 
our posterior In particular we work with an adaptive Metropolis algorithm utilizing a mixture proposal kernel 
known to satisfy these two ergodicity conditions for unbounded state spaces and general classes of target posterior 
distribution, see Roberts and Rosenthal [25 1 for details. 



5.0.2. Euclidean and Riemann- Manifold Adaptive Metropolis within Data- Augmented Copula PIC Models 

This section presents the specific details of the Adaptive Metropolis algorithm that we combine with Data- 
Augmentation to obtain an MCMC sampler for the Data Augmented Mixture Copula PIC Model proposed. This 
involves specifying the details of the proposal distribution in the AdMCMC algorithm which samples a new proposed 
update vector T* and matrix S* from an existing Markov chain state T with 

T = '4', Sq. j, tQ. j, p, Xi^^^^, . . . , Xi^^y^^ . . . , -X'j jjjj^ 

and matrix S. At the j-th iteration of the Markov chain we have existing state Y^^~^) and which is used to 

construct the proposal distribution q (T'-'"^^ , Y*) q ^E^^~^\ E*^ . The choices we make for the two proposals will 
involve a novel development of a new adaptive proposal for positive definite matrices, required for the covariance 
matrix S should we choose not to specify it as diagonal. 
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Euclidean Space Adaptive Metropolis for Static Parameters: 

Wc first detail the proposal for updating T using a mixture of multivariate Gaussian distributions as specified for 
an Adaptive Metropolis algorithm which involves sampling from the proposal 

q (t(*-i), .) = w^N- (t; Y(*-i), i^Cov (^{y«}^^^^^_ + (1 - w^)M (t; t(*-i) , Hi!/, , j , (5.1) 

where we define the sample covariance for Markov chain past history by Cov ^{T^^^ j^^^^.^^ .^^^ and we note the 
following recursive evaluation, which significantly aids in algorithmic computational cost reduction 



\ I. )0<j<t-lj 



(5.2) 

The theoretical motivation for the recommended choices of scale factors 2.38, 0.1 and dimension d are provided in 
Rosenthal et al. 26[. 



Riemannian Manifold Adaptive Metropolis for Covariance Matrices: 

Next we develop a novel proposal distribution for the sampling of the covariance matrix S S Sym^ (d) in an adaptive 
MCMC proposal, restricted to the Riemann manifold of symmetric, postive definite {d x d) matrices, denoted by 
the space Sym^(d). 



its distrihu- 



Remarks 5.1. First, we note two properties of the marginal posterior p E 

\ I. 0<i<J y 

tion is restricted to the Riemann-manifold of symmetric positive definite matrices, but in general will not be Inverse- 

Wishart; second, the Markov chain samples drawn from this marginal distribution at iteration t, < S'^-* > , are 

I J 0<s<t 

not independent. The consequence of this is that we cannot simply apply the property of closure under convolution 
of independent Wishart distributed random matrices to find a suitable proposal. 

Therefore, we will adopt a strategy to perform adaptive moment matching of a distribution with support Sym^(c?). 
We detail one possibility involving an inverse Wishart distribution fitted to the sample mean of the marginal 
posterior for the covariance. We note that future work could also consider specifying a distribution on the superset 
of the Riemannian manifold of symmetric positive definite matrices, given by the Riemannian manifold of symmetric 
matrices Sym^((i) C Sym((i). 

Adaptive Metropolis inverse Wishart Mixture: We note that one way to achieve this is a mixture of inverse 
Wishart distributions given by 

(^E;Af'^^' ({^^'^}o<.<*-i) ^ p) + - {^-^ p) ■ (5-3) 

Here, the adaptive proposal mixture component is specified through fixing the degrees of freedom p and then 

selecting A?''"^ with respect to the sample average of the covariance matrices S'^'*^ > which are samples 

L J o<s<t-i 

from the matrix-variate marginal posterior in the Markov chain, thereby adapting the proposal to the Markov 
chain history. To perform the moment matching (Equation (j5.4l) ). we note that wc need to ensure that the sample 
average considered is restricted to the Riemann-manifold of positive definite matrices. 

0<s 



^^"1{^'^U<.-J-^ ' (.-d.^(i:)-l). (5.4) 
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This is satisfied through the choice of the estimator 

S (5.5) 

s=l 

To see this we observe that since we only form positive linear combinations of matrices on this manifold, with a 
scaling, such linear combinations will always remain on the manifold Sym^(cZ). 



6. Real Data Analysis 



To illustrate the proposed models and compare with existing models and estimation methods in the actuarial 
literature we consider, as in Merz and Wiithrich the example presented in Dahms jij and Dahms et al. Q 
(Tables 10 and 11). As in the second analysis framework in Merz and Wiithrich [l^, we treat the claim development 
factors, the likelihood dependence parameters and the hyperparameters on the claim development factor priors as 
parameters which wc incorporate into the posterior inference. 

We present two sets of results, the first studies the performance of the adaptive Markov chain Monte Carlo 
algorithms developed for the estimation and inference of the posterior distributions for the PIC-Copula models 
for Gaussian Copula (Models III) and the Data-Augmented-Mixture-Copula PIC (Models IV). The second stage 
of results assesses the estimation of predictive distributions and dependence features of the PIC claims reserving 
models compared to the independent PIC Model, the payment only model and the incurred only models. In 
particular, we focuss analysis on the data sets studied in Merz and Wiithrich [l^ for comparison of the infiuence 
of dependence features in PIC models versus independence assumptions when performing PIC claims reserving. 

Convergence Analysis: In all the Markov chain Monte Carlo simulations, for each model (payment, payment- 
incurred Gaussian copula Model III; and Data- Augmented hierarchical Archemdean mixture copula Model IV), 
we carried out convergence diagnostics. This included the Gelman- Rubin R-statistics (all less than 1.5), the ACF 
plots for each parameter were checked to ensure all parameters had AGF's which were less than 10% by lag 20. 
Then the first 20% of samples were discarded as burnin and the remaining samples were used in inference results 
presented below. 



6.1. Results: Euclidean and Riemann-Manifold Adaptive Metropolis for hierarchical Bayesian Copula PIC Models 

In the simualtion results, we consider a block Gibbs sampler with the following three stages performed at each 
iteration of the adaptive Metropolis-within-Gibbs sampler for the PIC Model III and Model IV: 



Stage 1: Perform exact sampling of the development factors and their hyperparameters under the conjugacy 
results developed. 



Stage 2: Perform Euclidean space Adaptive Metropolis updates of the Augmented Data variables using 
proposal in Equation (|5.ip . 
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stage 3: (Gaussian Copula Model III) - Perforin Riemannian space Adaptive Metropolis updates of the 
covariancc matrix in the Gaussian copula. Note, we consider the constrained specifications presented in 
the "Dependent Lag Years" model specification in Section 13.21 Equation p.23p . Under this hierarchical 
Bayesian model, the joint covariancc between all observed payment and incurred loss data under the dependent 
development years assumption, satisfies a telescoping diagonal block size form covariancc matrix structure. 
Hence, the sampling of this structure can be performed blockwise on each covariance sub-block; 
(Mixture Clayton-Gumbel Copula Model IV) - Perform Euclidean space Adaptive Metropolis updates of the 
mixture copula parameters. 



6.1.1. Hierarchical Bayesian Gaussian Copula (telescoping block covariance) PIC (Model III) 

This section presents the estimation results for the Gaussian Copula based PIC models (Model III) on the real 
data. Figure [5] summarizes the dependence structure by a heatmap for the posterior distribution of the Gaussian 
copula covariance matrix. As mentioned in the introduction, the telescoping block covariance refers to the fact 
that the covariance structure is reducing in rank by 1 on each diagonal block for the payment data and then the 
incurred data. This model has the joint covariance between all observed payment and incurred loss data under the 
assumption that the development years are dependent, satisfying a telescoping diagonal block size form covariance 
matrix structure. Summarising the information from such posterior samples for distributions of covariance matrices 
is non-trivial as discussed in Tokuda ct al. |28ll. where they develop a four layer approach. Our article adopts aspects 
of the ideas proposed in Tokuda ct al. [23] to interpret the features of the posterior distribution samples for the 
dependence structures. 

The posterior mean for estimated PIC covariance structure is obtained by using Monte Carlo samples from the 
Ricmann-Manifold Adaptive Metropolis sampler and given by the estimator, 

s ( / J \ /J 



E 



E|P,/ 

(^) 




p 



e^n , (6.1) 



.i=0 



where | © {®i=o ^o) | i'' the s-th sample of the J( J — 1) x J( J — 1) covariance matrix. The 

estimated posterior mean covariance matrix is reported in a heatmap for the correlation matrix in Figure [5] In 
addition, we present examples based on posterior mean covariance for covariance sub-blocks p {Y.^\P, /) and then 
for p (E4IP, /), where e SP^{&) and S4 G S'P+(5), again converted to heatmaps of the correlation. We see 
that although the priors selected for the dependence features in Model III in all cases favoured independence, since 
the scale matrices were all diagonal i.e. Al' = Ig and A| =15, the resulting summaries of the marginal posteriors 
of the covarianccs clearly indicate non-trivial dependence patterns in the development years within the payments 
data and the incurred loss data. This is observed throughout each sub-block covariance matrix. 

Table [5] provides a second summary of the posterior for the covariance matrix which further demonstrates features 
of the dependence properties in the payment and incurred data per accident year and involves the estimates of the 
largest eigenvalue of each block diagonal matrix for the payment and incurred data as summary statistics. These 
estimates are given by 

a|'^ = argmax (det(E|'^ - Al) = o) . (6.2) 
The largest eigenvalue provides information on the posterior distribution of the magnitude of the first principal 
component of each development year, decomposed by accident year. That is, we can quantify in the PIC model, by 
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Figure 2: Top panel: Heatmap of the posterior distribution for the Gaussian copula covariance matrix (100 X 100), summarised by 
the heat map for the mean of correlation structure using samples from the Riemannian Manifold Adaptive Metropolis sampler under 
restriction to a telescoping diagonal block form. Bottom Left Panel: Heatmap for the posterior distribution sub-block covariance 
matrices Sq' and converted to correlation matrices. Bottom Right Panel: Heatmap for the posterior distribution sub-block 
covariance matrices Sq and S( converted to correlation matrices. The color key is given at the top left. 

accident year, the proportion of residual variation in the log payments for accident year i currently unexplained by 
the development factors $o:,/-ii which were jointly estimated in the PIC model and assumed constant accross each 
accident year (i.e. constant per development year) for parsimony. We can also repeat this for the incurred loss data. 
Suppose that a principal component analysis is performed, decomposing the variation in the payment and incurred 
data for each accident year i with respect to the variation unexplained by the development factors in the PIC model. 
Then, up to proportionality, the distribution of the eigenvalues corresponds to the proportion of contribution from 
the leading eigenvector (principal component). When this is coupled with the fact that we can also easily obtain 
samples from the marginal posterior distribution of the leading eigenvector of the covariance matrix for the i-th 
accident year's payment of incurred loss data in the PIC model, then we get complete information per accident 
year on the ability of the development factors in the PIC model to explain variation in the observed loss data. 
Table [ni summarises the results for the average PCA weight (largest eigenvalue) and average posterior eigenvector. 

Tokuda et al. develops a framework which formalizes an approach to the summary of dependence structures. 
For the running example of results that we present for distributions p (E^'ji-', j) and p (EljP, j), under such an 
approach the third and fourth layers of summary are presented in Figure |31 This involves the presentation of 
contour maps of these marginal posteriors that are constructed using adaptive MCMC samples of these matrices. 
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In Figure |H the development factors for payment and incurred data marginal posterior distributions are presented 
along with the posteriors of the hypcrparameters for the Gaussian Copula based PIC models (Model III). Finally, 
we also compare the estimated posterior marginal distributions of the development factors for the payment and 
incurred loss triangles for the models: payment only model; the incurred only model; the Gaussian Copula (Model 

III) dependent model; the PIC [Full] independent model and the PIC [Partial] independent model of Merz and 
Wiithrich [l^- The results of this comparison include the posterior mean estimates of E ['I'ilP, /] and E [5'i|P, /], 
for all i G {0, 1, . . . , J} and the posterior quantilcs for left and right tails as measured by the fifth and ninety- 
fifth percentiles, given in Table [T] We note that the results in this section for the Gaussian copula models are 
obtained using the log ratio observational data and the restults for the Mixture Archimedian copula model are more 
conveniently obtained using the log observations (not ratio data). 

It is also worth noting other approaches that can be adopted in the case of the Gaussian copula model. One 
could also included a data-augmentation stage in the analysis as was utilised in the Mixture Archimedian copula 
example. In addition, the covariance matrices could have been specified under different structures with more or 
less parsimony. The examples utilised in this section were those which provided a reasonable tradc-off between 
parsimonious model specification, while allowing a meaningful decomposition of the results. 

The results of the comparison between the Gaussian copula PIC model and the independent PIC model illustrated 
that whilst the posterior marginal mean development factor estimates are not affected by the dependence feature 
included, the marginal posterior shape is affected. This is reflected by the comparison of the posterior confidence 
intervals for the Gaussian copula PIC model when compared to the payment or incurred individual models where 
there is a significant difference present in the shapes of the marginal posterior. It is expected that this will have 
implications for the estimation of reserves using these different will be quantified in the next section. 

6.1.2. Data- Augmented hierarchical Bayesian Mixture- Archimedian Copula PIC (Model IV) 

This section presents the estimation results for the mixture of Clayton and Gumbel Copula based PIC models 
(Model IV) on the real data are presented in this section. Figure [S] presents a summary of the mixture copula 
dependence structure obtained from posterior samples of the copula parameters under the hierarchical Bayesian 
model. The results in this section arc obtained using the log observational data, not ratio data. The figures 
summarise succinctly the estimated posterior dependence structure for the hierarchical Bayesian mixture Copula 
model, through plots of the dependence structure as captured by the estimatd mixture copula distribution, the 
scatter plots of copula parameter for the lower tail and rank correlation (Kendall's tau) and the upper tail copula 
parameter versus rank correlation. These results clearly demonstrate posterior evidence for non-trivial tail depen- 
dence features in the payment and incurred data, as well as potential for asymmetry in the upper and lower tail 
dependence. Note, uniformative prior choices were made on the copula parameters with uniform priors over [0, 50] 
and [1, 50] respectively, indicating these estimated copula parameters are data driven results. 

Figure [H] presents the development factors for payment and incurred data marginal posterior distributions along 
with the hypcrparameter marginal posteriors for the Data- Augmented Mixture Copula based PIC models (Model 

IV) . 
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Figure 3: Hcatmaps for the block diagonal covariance matrices T,^ (2x2 sub-plot 1) and E| (2 X 2 sub-plot 2). These are obtained 
using samples from the Riemannian Manifold Adaptive Metropolis sampler. Samples from the Posterior distribution of the telescoping 
diagonal block size form covariance matrix structures of the Gaussian copula under the hierarchical Bayesian model which has the 
joint covariance between all observed payment and incurred loss data under the dependent development years. Each set of 4 X 4 
panels, starting from the top, summarizes the posterior distributions for the covariance matrices for s G {P, 1} according to: Top Left 
Panel: contour map of posterior samples log [S|] j vs log [S|] ^ j. . Top Right Panel: contour map of posterior samples log [S|] ^ ^ 
vs g. Bottom Left Panel: kernel density estimator of the posterior distribution of the trace of the covariance matrix using 

samples {logtr Bottom Right Panel: scatter plot of posterior samples of the first, second and third largest eigenvalues scaled 

by total of the eigen valusc - (PCA weights - for linear combinations of the development factors when explaining variation in observed 
payment and incurred data for a given accident year). 
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Figure 4: Boxplot summaries of the marginal posterior distributions obtained using samples from the Riemannian Manifold Adaptive 
Metropolis sampler. Samples from the Posterior distribution under a telescoping diagonal block size form covariance matrix structures 
of the Gaussian copula under the hierarchical Bayesian model which has the joint covariance between all observed payment and incurred 
loss data under the dependent development years. Top Left Panel: box plots of marginal posterior distributions for p{^i\P, I). Top 
Right Panel: box plots of marginal posterior distributions for p('J'i|P,7). Bottom Left Panel: box plots of marginal posterior 
distributions for p {si\P, I). Bottom Right Panel: box plots of marginal posterior distributions for p {ti\P, I). 
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Figure 5: Copula Dependence Parameter Posterior distributions estimated under tlie Data- Augmented Mixture Copula PIC Model 
IV. A mixture of Archimedean copula models is considered, with Clayton and Gumbel copula choices, allow for possible asymmetry 
in the tail dependence over development years. We chose uniformativc uniform priors U[0, 20] for the copula parameters. Top Left 
Panel: Contour map of posterior estimated mixture copula dependence distribution between development years over paid and incurred 
loss data, with homogeneous dependence assumptions over accident years (estimated from posterior mean of p^'^'^^ and pQ^^^^ . 
Top Right Panel: Surface plot of posterior estimated mixture copula dependence distribution between development years over paid 
and incurred loss data, with homogeneous dependence assumptions over accident years (estimated from posterior mean of p^^^^^ 
and pQ^^^^ . Bottom Left Panel: Scatter plot of posterior samples used to estimate Kendall's tau rank correlation versus copula 
parameter for the Clayton mixture component. Bottom Right Panel: Scatter plot of posterior samples used to estimate Kendall's 
tau rank correlation versus copula parameter for the Gumbel mixture component. 
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Figure 6: Boxplots of the marginal posterior distributions of the development factors and hyperparameters. Top Left Panel: box plots 
of marginal posterior distributions for p {^i\P, I). Top Right Panel: box plots of marginal posterior distributions for p(^i|J-',7). 
Bottom Left Panel: box plots of marginal posterior distributions for p {si\P, I). Bottom Right Panel: box plots of marginal posterior 
distributions for p {ti\P, I). 
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7. Comparison of PIC reserving with Gaussian Copula PIC and Mixture Archimedian Copula PIC 
Models 



This section discuss the effect of modening the dependence structures on the reserving estimates. First we note two 
important details in calculating the reserves. We need to be able to draw samples from the predictive distributions 
for the payment and incurred data given below, for each accident year i, using 

p(P,,j|P,7)= J p{P,j\P,^i.,j_,,@)p{@\P,I) d@ and pihj\P,I)= J pihj\h^i.,j^„@)pi@\P,I)d&. 

In general it is not possible to solve these integrals analytically. Howerver, for the Gaussian copula models developed 
in this paper, under the results in Lemma [3.101 one adopt the results of Merz and Wiithrich [l^ [Theorem 2.4] to 
obtain analytic Gaussian predictive distributions. Alternatively, the predicitive distributions can be estimated as 
described in Peters et al. Section 3.3]. Although the results in Table [T] demonstrate that the incorporation 
of the dependence structures does not significantly alter the posterior mean of the development factors for the 
payment and incurred loss data, it is clearly possible for the predictive distribution to be altered, since the shape 
of the posterior distribution is altered by the dependence features. Second, regarding the hierarchical mixture 
Archimedian copula model, it does not admit an analytic solution for the predictive distribution. This does not 
matter if a data augmentation stage is set up in the joint posterior distribution to sample cumulative payments, 
since then we can use the MCMC sampler output for the ultimate cumulative payment and incurred losses in each 
accident year. 

Finally, we also note that a simple Monte Carlo based approximation for the ultimate claim can be constructed. 
Take the samples from the MCMC output for the PIC model of interest (sampled from the complete PIC model with 
dependence features present) and then utilise these samples to construct a Laplace approximation to the predictive 
observation distribution for example p {Pi,j\Pi,i:j-i, &) which involves a normal approximation around the MAP or 
locally around each Monte carlo sample for the development factors, with precision given by the sampled observation 
covariance structure. Though this is not required, as we have shown for the Gaussian copula models independence 
models, it may be useful for alternative copula based models with simple data-augmentation approaches. In 
addition a second alternative would be to utilise in the predictive distribution the marginal distributions. 

Figure [7] presents the log posterior predictive distribution for the ultimate total claim given by the predictive 
distribution for the log of the cumulative payment over each accident year X]f=o '^'-'^ ^^^^ Bayesain models 
which incorporate priors on observation error, development factors and hyperpriors for precision of the development 
factors. We see that all three models are in good agreement with each other with the dependence parameters 
affecting the variance and tail behaviour of the distributions. 

Next we consider the distributions of the outstanding loss liabilities estimated using the S samples from the 
MCMC obtained for the posterior PIC model. We denoted these by random variables {i?(P, J)^*-' where 
— Pi, J — Pi^j-i and depending on whether payment, incurred, or both data is present we denoted by 
i?(P)(^\ Ril)'-"'^ and R{P,I)'-^^ respectively. Figure [HI presents the MCMC estimated claims reserve marginal 
posterior predictive distributions for each accident year per model developed. We compared our results to those 
obtained in Merz and Wiithrich and find good agreement between the mean reserve per accident year and 
each proposed model. In addition, we note the possible differences between the distributions can be attributed to 
the utilisation of the full versus partial hierarchical Bayesian models in this paper and the different dependence 
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PIC Gaussian Copula (Full) - \- - 



-]1 



PIC Independent (Full) - 



Payment Only (Full) - 



PIC Mixture Copula (Full) - 



16 16.5 17 

Log Ultimate Total Payments 



PIC Gaussian Copula (Full) 
PIC Independent (Full) 
Payment Only (Full) 
PIC Mix Copula (Full) 




Figure 7; Boxplots of the predictive distributions obtained from the MCMC samples. Ultimate Bayesian predictive distributions for 
log payment data from the payment only predictive distribution, the Full Independent PIC model, and the hierarchical PIC Mixture 
Copula model via Data Augmentation predictive distribution. Left Panel: Posterior predictive distribution box plots from samples. 
Right Panel: Kernel density estimates of the predictive distributions. 



structures considered. Additionally, we note that further analysis on comparisons to existing models in the literature 
can be obained for the models of Mack [l^ , Dahms Q and Quarg and Mack 

and Wiithrich Q [Table 4] and in the spreadsheet provided by Professor Mario Wuethrich at URI^. 



2J] for this data analysis in Merz 



8. Conclusions 

This paper extends the class of PIC models to combine the two different channels of information as proposed in 
Merz and Wiithrich 



16| by introducing several novel statistical models for the dependence features present within 
and between the payment and incurred loss data. This allows us to obtain a unified ultimate loss prediction which 
incorporates the potential for general dependence features. To achieve this we developed full hierarchical Bayesian 
models which incorporate several different potential forms of dependence, including generalized covariance matrix 
structure priors based on inverse Wishart distributions and conditional Bayesian conjugacy in the PIC independent 
log-normal model. This forms a general class of Gaussian copula models which extends the approach of Happ and 
Wuthrich [ll|. 

Second, we develop a class of hierarchical mixture Archimedian copula models to capture potential for tail depen- 
dence in the payment and incurred loss data, again developing and demonstrating how to appropriately construct 
a full Bayesian model incorporating hyperpriors. In this regard, we also develop a class of models in which data- 
augmentation is incorporated to both overcome challenging marginal likelihood evaluations required for the MCMC 
methodology to sample from the PIC Bayesian models. This had the additional feature that it also allowed for 
joint Bayesian inference of the reserves as part of the posterior inference. 



^ URL fhttp : //www ■ math . ethz ■ ch/- gueth/clalms_reserving3 . html | 
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Figure 8: Boxplots of log ultimate Baycsian predictive reserve distributions for payment data per accident year, compared to (Partial) 
PIC Independent posterior mean estimates from Merz and Wiithrich [jj (kargc unfilled black circles). Top Row: the (Full) hierarchical 
PIC Mixture Copula model via Data Augmentation; Second Row from Top: the (Full) hierarchical PIC Gaussian Copula model; 
Third Row from Top: the (Full) Independent PIC model; Bottom Row: the (Full) payment Only model. 



Finally, to perform inference on these approaches we developed an adaptive Markov chain Monte Carlo sampling 
methodology incorporating novel adaptive Ricmann-manifold proposals restricted to manifold spaces (postive def- 
inite symmetric matrices) to sample efficiently the covariance matrices in the posterior marginal for the Gaussian 
copula dependence. We make these advanced MCMC accessible to the actuarial audience to address challenging 
Bayesian inference problems in Claims Reserving modelling. 

The consequence of these models for actuaries is that a new extended suite of flexible dependence structures have 
been incorporated into the recently proposed PIC models. These can now be extended and compared to existing 
chain ladder methods. We perform an analysis on real payment and incurred loss data discussed in Merz and 
Wiithrich [l^ and compare our models with the analysis for the independent PIC model (partial) and the (full) 
Bayesian PIC model as well as several different dependent models and the payment only model. Furthermore, we 
jrovide reference on further comparisons to the alternative models of Mack 14 1, Dahms ^ and Quarg and Mack 



2J] for this data. 
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Posterior Covariance Matrix for payments and incurred Loss Gaussian Copula 
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Posterior Marginal Distributions for Development Factors 
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Table 1: NOTE: (Pull) corresponds to PIC models with results for the FULL hierarchical Bayesian PIC model with priors on development factors, observation variances and hyperpriors on 
precisions on development factors. The PIC Independent (Partial) of Merz and Wiithrich are the Bayesian posterior results in which a-m and t„ are assumed known. In addition, the PIC 
Mixture Copula model has posterior development factors on the scale of log cumulative payment data (not ratio data), so the reported posterior mean development factors are for the cumulative 
payment marginal posterior means (log scale). 



A. Appendix A 



Lemma A.l (Properties of Matrix- Variate Gaussian Distribution). A p x n random matrix X is said to 
have a matrix variate Gaussian distribution with p x n mean matrix M and covariance matrix E (g) ^' where S 
and \E' are in the spaces of symmetric positive definite matrices given hy Yi £ SD"*" {MP) and \1/ € SD"*" (K") if the 
pnx 1 dimensional random vector Vec {X') has a multivariate normal distribution Vec {X') ~ N (Vec{M'), S ® 5'). 
Furthermore, if X is distributed according to matrix-variate Gaussian distribution X ^ Np_„ (M, E Cg) 4') then the 
density is given by 

exp {-^tr [S-i {X - M)' -^-^ {X - M)]) 



(27r) 



np/2 |j]|"/2 



(A.l) 



In addition the following properties are satisfied for a matrix-variate Guassian (see ^] Ghapter 2): 



1. If X 

2. IfX 



Np^n {M, S (g) «-), then X' - 7V„,p (A'/', 

Np^n (Af, S ® *),, and partition X, M, E, and * as 



X 



Xir 
X2r 



and X = [Xic X2c] 



(A.2) 



with Xir the (m x n) sub-matrix, X2r the [p — m, x n) sub-matrix, Xic the {p x t) sub-matrix and X2c the 
ipxn — t) sub-matrix. With analogous partitions of the mean matrix Mir, M2r, Afic and M2c and covariance 
matrices 



S21 E22 



and * = 



*11 1'12 
1'21 1'22 



(A.3) 



with Ell the (m x m) sub-matrix, S12 the (m x p — m) sub-matrix, E22 the (p — m x p — m) sub-matrix, '^n 
the {t X t) sub-matrix, '^22 the (n — t x n — t) sub-matrix. Then the following properties are true 



Xir ~ A^m,n (Afi^, En ® ^) and Xi^ - Np^t (Afic, E ® *ii) 
X2^|Xi^ - A^p-„,,„ (Af2^ + EaiEj;/ (ATi^ - Mir) , S22.1 ® 

X2c\Xic - Np^r.-t {M2c + {Xi, - Ml,) *ri'*12, ^ ® «'22.l) 



(A.4) 



with E22.: 



E2iEii^Ei2 anrf ^-22.1 = ^'22 - *2i*ii^*i2. 



Lemma A.2 (Properties of Matrix Variate Wishart Distributions). Suppose that the {p x p) positive def- 
inite matrix E has an inverse Wishart distribution, with positive definite (p x p) scale matrix A, degrees of freedom 
parameter k > p — 1, and density 



/(S|A,fc) 



|j^|-(fc+p+l)/2 g-tracc(AS-i)/2 



(A.5) 



2'=p/2rp(fc/2) 

where Tp{-) is the multivariate gamma function. The mean and mode of this distribution are given respectively by 



E[E|A,fc] = 



k — p — 1 



A, and m (E) = 



1 



k+p+1 



-A. 



(A.6) 



Furthermore, the following marginal and conditional properties of the inverse Wishart distribution are relevant. 
Gonsider a partition of the matrices A and as 



A 



All A12 
A21 A22 



Sii E12 
S21 E22 



(A.7) 



with Aij and Ej^- denoting pi x pj matrices, then the following properties are satisfied (see Gupta and Nagar 
f^][Ghapter S, Section 3.4]): 

1. The random sub-matrix En is independent o/Ej~^^Ei2 

2. The marginal distribution of any sub matrix on the diagonal of the matrix E is distributed as inverse Wishart. 
For example, the sub random matrix En is as inverse Wishart with En ^ IW (An, k — P2); 

3. The marginal distribution of sub random matrix E221 is inverse Wishart E221 ^ IW (A22.1, fc). 



In Lemma IA.3I below we present details for the matrix-variate Inverse Wishart distribution, see 9| [Chapter 3.4, 
Definition 3.4.1 and Theorem 3.4.1] 
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Lemma A. 3 (Properties of Matrix- Variate Inverse Wishart Distribution). A random p x p matrix V = 
is distributed as Inverse Wishart, with degrees of freedom m and p x p parameter matrix '5, denoted V ~ 
IWpljn,'^) with density 

fim, m) = r fl/2fm-J- nily|i/^"' "^^ , F > 0, * > 0, m > 2p. (A.8) 



B. Appendix B 



The family of Archimedean copula models has the following useful properties presented in Lemma IB. II 

Lemma B.l. Let C he an Archimedean copula with generator ip. Then according to Nelsen Lemma 4.1-2 and 
Theorem 4-L5], the following properties hold: 

1. C is an Archimedean copula if it can be reprsented by 

C{u, v) = (^["^l {if{u) + if{v)) 

where ip is the generator of this copula and is a continous, strictly decreasing function from [0, 1] to [0, oo] 
such that = and ip^~^^ is the pseudo inverse of ip. 

2. C is symmetric, C{u, v) — C{v, u) W{u, v) G [0, 1] x [0, 1] 

3. C is associative, C(C(u, v),w) = C{u, C{v, w)) V(m, v, w) e [0, 1]'^. 

4. If c > is any constant, then cip is a generator of C 

5. According to Denuit et al. 0, Definition 4.7.6], the extension of the Archimedean copula family to n- 
dimensions is achieved by considering the strictly monotone generator function ip such that p : (0, 1] — > 
with p{l) = 0, then the resulting Archimedean copula can be expressed as, 



C(ui,U2, . . . , 




The members of the Archimedean copula family utilised in this manuscript are given below in Lemma IB. 21 



Lemma B.2. From the results in Nelsen fl8i . Section 4.3, Table 4.1] the distribution and density functions of the 
Clayton copula are given respectively as: 



C'^{ui, Un) = 1 - n + 



(B.l) 



(B.2) 



where p'^ € [0, oo) is the dependence parameter. The Clayton copula does not have upper tail dependence. Its lower 
tail dependence is \l = l^^l f . The distribution function of the Cumbel copula is 



C^{ui, ...,Ud) = exp 



.1=1 



(B.3) 



where p'~' G [l,oo) is the dependence parameter. The Gumbel copula does not have lower tail dependence. The 
upper tail dependence of the Gumbel copula is Xjj = 2 — 2^^^ . The distribution function of the Frank copula is 



C^{ui, ...,Un) 



p y {eP'^ - !)"-! J 



(B.4) 



where p^ G M/{0} is the dependence parameter. The Frank copula does not have upper or lower tail dependence. 
We note that the density functions for Cumbel and Frank does not admit simple recursive expressions in terms of 
their density functions, but they can be obtained via partial differentiation 

gn 

c{ui,...,Un) = C(ui, ...,M„). (B.5) 

OUl, ...OUn 
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C. Appendix C 



Proof The proof of Lemma \4.1\ requires one to demonstrate that the resulting distribution function 

J[0,tli]x[0,tl2]x---X[0,M„] 

^ Wi Ci{xi,X2, ■ ■ ■ ,Xn) dxi.,n 



[0/ui]x[0,ti2]x---x[0,ti„] 



= ^ WiCi {ui,U2, . . . ,u„) 

satisfies the two conditions of a n-variate copula distribution given in [Definition 2.10.6] of Nelsen [id /. The first 
of these conditions requires that for every u — {ui,U2, . . . , Un) G [0, 1]", one can show that C (u) = if at least 
one coordinate of u is 0. Clearly since we have shown that C (u) = X^j^i"" '"^iC'i {iJ') given each member 
Ci {ui,U2, . . . ,Un) G yl" is define to be in the family of Archimedean copulas each of which therefore satisfies 
this condition for all such points u, then it is trivial to sec that the probability weighted sum of such points also 
satisfies this first condition. Secondly one must show that for every a and b in [0, 1]", such that a < b (i.e. 
ai < bi Vi £ {1, 2, . . . ,n}) the following condition on the volume for copula C is satisfied, V(j{[a, b]) > 0. As in 
Nelsen we adopt the notation for the n-box, [a, b], representing [ai, bi] x [02,62! x • • • x [o-m Ki] 8,nd we define 
the n-box volume for copula distribution C by [Definition 2.10.1, p.43] of Nelsen \lH] giving 

Vc{[a,b]) = J2sgn{c)C{c) 

where the domain DomC of the mixture copula C satisfies [a, b] C DomC. In addition we note that this sum is un- 
derstood to be taken over all vertices c of n-box [a, b] and sgn{c) = 1 if Ck ~ o-k for an even number of k 's or sgn{c) = 
— 1 if Ck = dk for an odd number of k 's. Equivalently, we consider A^^^C (t) = C* (ti, ^2, • • • , tk-i, bk, tk+i, ■ ■ ■ ,tn) ^ 
C {ti,t2, ■ . ■ , tk-i, flfc, ifc+i, . ■ . ,tn). In the case of the mixture copula, we can expand the volume of the n-box [a, b] 
as follows 

m ra 

Vci[a, b])=J2 sgnic)C (c) = ^ ^ w,sgnic)Q (c) = ^ ^ m.^c, ([a, b]) 

i=l 1=1 

hence we sec that since each component Ci (ui, U2, . . . , w„) is a member of the set of Archimedean copula distribu- 
tions A'^, therefore for each component we have that Vci{[a, b]) > for all i £ {1, 2, . . . , m}. 

□ 
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